Failure of mean-field approach in out-of-equilibrium Anderson model 



O 

o 

O 

Q 



i 

C 

o 
o 



> 

OS 
(N 

o 
o 



X 



B. Horvath 1 , B. Lazarovits 1,2 , 0. Sauret 1 , G. Zarand 1 
1 Theoretical Physics Department, Institute of Physics, 
Budapest University of Technology and Economics, Budafoki ut 8, H-1521 Hungary 
2 Research Institute for Solid State Physics and Optics of the Hungarian Academy of Sciences, 
Konkoly-Thege M. ut 29-33., H-1121 Budapest, Hungary 

To explore the limitations of the mean field approximation, frequently used in ab initio molecular 
electronics calculations, we study an out-of-equilibrium Anderson impurity model in a scattering 
formalism. We find regions in the parameter space where both magnetic and non-magnetic solutions 
are stable. We also observe a hysteresis in the non-equilibrium magnetization and current as a 
function of the applied bias voltage. The mean field method also predicts incorrectly local moment 
formation for large biases and a spin polarized current, and unphysical kinks appear in various 
physical quantities. The mean field approximation thus fails in every region where it predicts local 
moment formation. 
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The Anderson impurity model^ (AIM) has been the 
subject of great theoretical and experimental interest in 
the past decades (for a review see Ref. d). There is a 
number of experimental systems including quantum dots, 
or single atoms and molecules contacted by leads, which 
provide experimental realizations of various versions of 
the AIM under out-of-equilibrium conditions. These sys- 
tems are not just prototypes of out of equilibrium sys- 
tems but a theoretical understanding of them would be 
crucial for future molecular electronics and mesoscopic 
applications. 

Anderson constructed his famous model in Ref. [l] to de- 
scribe local moment formation and solved it within the 
mean-field (MF) approximation. Within this approxi- 
mation, he found a phase transition to a state where 
magnetic moments are formed. Further work revealed 
that, in reality, quantum fluctuations of this local mo- 
ment lead to the formation of a Kondo-singlet between 
the impurity and conduction electrons 2 - at low tempera- 
ture, where the impurity spin is thus completely screened. 
The spontaneous symmetry breaking found by Anderson 
is thus an artifact of the mean field approximation. Nev- 
ertheless, the MF treatment indicates clearly the regions 
of strong correlations, and it can also serve as a starting 
point for accurate approximations as in the local moment 
approach- (LMA) or interpolative perturbation theory 4 ^ 
(IPT). The latter approach can easily be generalized to 
non-equilibrium situations^ using Keldysh formalism 8 . 

In lack of more accurate methods, the mean field ap- 
proximation is also used in molecular electronics calcula- 
tions, where LDA or eventually Hartree-Fock equations 
are solved in a scattering state or Keldysh approach to 
describe moment formation^. However, it is not clear at 
all, how reliable these approximations are. The purpose 
of this paper is to shed some light on the weaknesses of 
the non-equilibrium mean field approach on the simplest 
possible test case, the out of equilibrium Anderson model, 
and to show where usual ab initio calculations should fail. 
Our conclusion is that the mean field approach fails qual- 
itatively and quantitatively essentially everywhere where 



it predicts local moment formation. Our study, which is 
based on the scattering state formalism, is complemen- 
tary to the recent work of Komnik and GogoliniS, who 
used a Green's function formalism to study the mean field 
equations of the non-equilibrium Anderson model. As we 
shall see, in the strongly correlated regions several arti- 
facts emerge such as non-equilibrium driven spontaneous 
symmetry breaking as well as multiple stable solutions 
which lead to the appearance of hysteresis. These in- 
stabilities are probably also parts of the reasons, why 
non-equilibrium IPT suffers from all kinds of instabili- 
ties. These instabilities were avoided in previous works 
by applying spin-independent approximations^ 1 - or using 
an interpolative self-energy 4 or both 12 . 

The non-equilibrium AIM Hamiltonian consists of four 
parts. The first part describes a single impurity level with 
energy Ed and an on-site Coulomb interaction (U) 
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where d\, and d a are the creation and annihilation op- 
erators of the impurity electrons corresponding to spin 
state a and n a = The second and third terms de- 

scribe the left (L) and right (R) leads which we model by 
tight-binding chains, 



H a = (— 2t cos k 



Here a £ (L, R), 



and 



Ckaa are the creation and 



annihilation operators of a conduction electron of wave 
number k € {0, 7r} in lead a, t is the hopping along the 
leads, and fi a is the chemical potential of the left or the 
right lead. The chemical potentials of the two leads are 
different due to a finite bias voltage leading to a non- 
equilibrium situation. The fourth part, H t , describes the 
tunneling between the leads and the impurity 



H t = VJ2 [4(V-C-la + V+c la ) + h.c 
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Here are the hybridization matrix elements between 
the impurity and the left and the right leads, and cj = ±i CT 
denote the conduction electron annihilation operators on 
the sites next to the impurity; sites along the left and 
right chains are labelled by / = {— oo, . . . , — 1} and I = 
{1, . . . , oo }, respectively. For the sake of simplicity, here 
we study a symmetrical situation, V- — V+ — V, but our 
conclusions are rather independent of this assumption. 

To study the Hamiltonian above we used a mean-field 
approximation and replaced the impurity term as 



H,, 



(4) 



The non-equilibrium expectation value of the occupation 
numbers, (n^), in Eq. (j4]) can be obtained by solving 
self-consistent equations discussed later. The expression 
£d+U (n-cr) can be regarded as an effective energy level of 
spin-CT electron. The hybridization between impurity and 
conduction electrons gives rise to a finite lifetime for im- 
purity states, reflected in the broadening of the effective 
impurity levels with a finite width, T = 2irV 2 po — V 2 /t 
where po is the density of states (DOS) of the conduction 
electrons at the Fermi-level of the half-filled leads. 

To evaluate the non-equilibrium expectation values, 
(n a ), we shall use a scattering formalism. The annihila- 
tion operator of the left-coming scattering state of energy 
e and spin a can be expressed as 
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Here ci a is the annihilation operator of the Zth site and 
a CT (e), Pa{s) and j a (e) are the reflection, transmission 
and dot coefficients, respectively, which we obtain by 
solving the corresponding Schrodinger equation. Since 
the energy is conserved in course of the scattering pro- 
cess, the wave numbers k and k' are connected by the 
dispersion relation e = — 2icosfc + pL = —2tcosk' + [ir. 
Waves coming from the right hand side can be con- 
structed in a similar way. The occupation numbers are 
then calculated from the dot coefficients, 7 CT (e) and Y a {e) 
corresponding to the states coming from the left and the 
right, respectively, 

/oo 
de(\ la (e)\ 2 fL(s) + WAz)\ 2 fn(e)) , (6) 
-oo 

where the DOS of leads is approximated by its value at 
Fermi-level, and f a {e) = f(e — p a ) is the Fermi function. 

In the large bandwidth approximation (t — > oo, T = 
V 2 ft finite) the occupation numbers become 
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FIG. 1: Magnetic (M), coexistence (C) and paramagnetic (P) 
regions as a function of p/U and T/U values for Sd/U = —1/2. 
The boundaries of the regions obtained by increasing and de- 
creasing bias are plotted with dashed and dotted lines, respec- 
tively. The exactly evaluated boundary curve is indicated by 
squares. Inset: Magnetization as a function of increasing and 
decreasing bias voltage for T/U = 0.05 (along the dotted line 
in the main figure) 



which simplifies further for zero temperature as 
1 ,_i ( £d + U{n- a ) - p a 



a£(L,R) 



r 



(8) 



These self-consistent equations can also be obtained us- 
ing the Keldysh formalisms^. Eqs. © or (jHJ) are solved 
iteratively by assuming that the applied bias is symmet- 
rical pl — m/2 and pn = —p/1. In these calculations the 
bias voltage is increased or decreased gradually and the 
local stability of the solutions is always checked. 

First, let us discuss the bias-dependence of the occupa- 
tion numbers at T = 0. In the inset of Fig.[T]the magneti- 
zation m = (n-f ) — (n j ) is plotted both for increasing and 
decreasing bias voltage at a fixed ratio of T/U = 0.05 
for the symmetric Anderson model (sd/U = —0.5). In 
equilibrium, at zero bias, the impurity possesses a finite 
magnetic moment for the chosen parameter set, while in 
case of high bias voltages (p > 0.5) the stable solution is 
paramagnetic. Between the two limiting cases a region 
appears, p c \ < p < p C 2, where both the magnetic and 
the non-magnetic solutions are stable. We shall refer to 
this region as a coexistence region. The existence of 
multiple stable solutions for the occupation numbers in 
this region is reflected in the hysteresis of the magneti- 
zation too. The sharp decay of the magnetization shown 
in the inset of Fig. [T] and the existence of a hysteresis 
between the critical fields indicate clearly a first order 
transition, predicted incorrectly by the MF solution. 

The parameter space can thus be divided into magnetic 
(M), paramagnetic (P) and coexistence (C) regions. In 
the paramagnetic regions a single stable solution exists 
only ((nj) = (n-i)), while in the magnetic one two sta- 
ble magnetic (corresponding to magnetizations ±m) and 
an unstable paramagnetic solution can be found. In the 
coexistence region two magnetic and one paramagnetic 
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FIG. 2: "Phase diagrams" as a function of n/U and T/U for 
T = and different £ d values, a) Ed/U = —0.35 6) —0.25, c) 
and d) 0.25. Notation of the regions and the lines are the 
same as Fig. Q] 



FIG. 3: "Phase diagrams" as a function of fi/U and T/U for 
e d /U = -0.5. The temperature a), T/U = 0.051 6), 0.101 c) 
0.151 and d) 0.201. Notation of the regions and the lines are 
the same as Fig. Q] 



stable solutions and two unstable magnetic solutions ex- 
ist. Therefore, these regions can be distinguished by the 
number of the solutions of Eq. ([8]). In the symmetric 
case when (n-j-) + (n^) = 1, the easiest way to construct a 
"phase diagram" is to sweep possible values of (rtj) — (n^) 
pairs, substitute them into Eq. ([5]) and count the number 
of solutions for different parameter sets. 

In the paramagnetic region one can exploit the fact 
that Eq. ([5]) has a non-magnetic solution for every possi- 
ble parameter set. Therefore one can substitute (rif) = 
(ni) into Eq. and search only for non-magnetic states. 
The region of stability for this solution can then be de- 
termined through a linear stability analysi o 10 ! 13 . 

Fig. [1] shows a typical a magnetic "phase diagram" 
as a function of fi/U and T/U for the symmetric non- 
equilibrium AIM at T = 0. At the "upper critical line" 
fi C 2 the magnetic solution becomes unstable. In the mag- 
netic case, the effective levels corresponding to different 
spins are not equally occupied and lie at different ener- 
gies. The magnetic solution becomes unstable when the 
value of the bias voltage reaches approximately the effec- 
tive energy of one of the two differently occupied levels. 

The critical line /i c i in Fig Q] marks, on the other- 
hand, the border of stable paramagnetic solutions. In 
the special case, ed — —U/2 the two spin occupa- 
tions are (n-\,\_) = 0.5 in the whole paramagnetic re- 
gion, and the magnetic boundary equation simplifies to 
^i 2 + 16r 2 — 8UT/n = 0. This analytical result is nicely re- 
produced by our numerical stability analysis (see Fig.[T]). 

The values of fi c i and /i C 2 are functions of T/U and 
Ed/U; for increasing T/U, the coexistence and magnetic 
regions disappear and only the non-magnetic solution 
survives. In Fig. [2] the "phase diagrams" can be seen 
as a function of [i/U and T/U at T = 0. Depending on 
the value of Ed we can distinguish four different regions: 
empty regime (ed > 0), mixed valence regime (ed ~ 0), 
local moment regime (-U < Ed < 0) and a doubly oc- 
cupied regime (ed < — U) which behaves similarly to the 
empty regime by particle-hole symmetry. 



In the local moment regime, shown in Fig. [T] and 
Fig. f2Ji, for Sd « —U/2, there is approximately one elec- 
tron on the impurity forming a local spin moment. In 
equilibrium, this finite magnetic moment is predicted on 
the dot below a critical value of T/U, and this moment 
is destroyed by a large enough bias voltage. The coex- 
istence region only appears in this local moment regime 
and vanishes above Ed ~ —17/4 (see Fig.[2>). 

In the empty regime (ed > 0), shown in Figs. &p 
and [2H, the equilibrium magnetization completely dis- 
appears. However, surprisingly, the MF solution pre- 
dicts the appearance of local moments for small T's and 
2ed < H < 2ed + U biases. This intriguing local mo- 
ment formation has a simple physical meaning: For large 
values of U/T and 2ed < H < 2ed + U, the bias volt- 
ages are large enough to inject an electron to the empty 
level, however, they are not large enough to overcome 
the Coulomb energy of injecting a second electron to the 
local level. Therefore, electrons pass through the dot one 
by one, and a fluctuating magnetic moment appears on 
the dot. Note that in this regime, the magnetization is 
induced exclusively by the finite bias voltage. 

The overall effect of the temperature in the applied 
MF-approximation is to destroy the magnetic moment on 
the impurity and drive the system to be paramagnetic. 
Fig. [3] shows the temperature dependence of the "phase 
diagram" in the symmetric case. The coexistence region 
gradually vanishes for increasing temperatures, while the 
magnetic and paramagnetic regions get larger. Increasing 
the temperatures further the magnetic region disappears 
too. 

The previous results are summarized in Fig. [4] for a 
fixed ratio T/U = 0.05 that is already sufficiently small 
to find the coexistence and magnetic regions in the local 
moment regime. The whole "phase diagram" is symmet- 
ric to ed/U = —0.5 due to the electron-hole symmetry. 
Note that the coexistence region only exists close to the 
electron-hole symmetry, ed/U ~ —0.5. For large asym- 
metries (double or zero equilibrium occupation) a non- 
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FIG. 4: Magnetic (M), coexistence (C) and paramagnetic (P) 
regions as a function of fi/U and Ed/U values, for T/U = 0.05. 



equilibrium magnetic solution appears, while in equilib- 
rium only the paramagnetic solution exists. This mag- 
netic region has also been observed although not dis- 
cussed in detail in Ref. [lfj. 

The mean-field solution also leads to the appear- 
ance of non-physical features in the transport proper- 
ties of the impurity. We calculated the transport prop- 
erties within the MF approximation using the Landauer- 
Biittiker formalism 1 ^, but similar results can be obtained 
applying the Keldysh formalism 1 ^. The current can be 
evaluated as 



la = 



- J Je(f L {e)-f R {e))\t a {e)\ 2 



(9) 



Here the transmission coefficient t a = /3 a \Jvy Jvk is nor- 
malized to the flux, with Vk and v^i the velocities of inci- 
dent and transmitted electrons with wave numbers k and 
k! . Applying the large bandwidth approximation again, 
the current can be written for finite temperatures as 



la = 



e 
2nh 



de- 



r 2 (Me)-f R (e)) 
(e-e d -U(n-a)f +T 2 



(10) 



with (ria) the non-equilibrium occupation numbers ob- 
tained from Eq. |(7J). This integral can be trivially evalu- 
ated at T = temperature. 

Fig. [5] shows the current as a function of bias for two 
different level positions in the strongly correlated regime. 
For Ed/U = —0.5 we find hysteresis in the current. For 
£d/U — 0.25, on the other hand, no hysteresis appears 
but the current shows a two-step behavior as a function 
of bias voltage. In this empty regime, the MF equations 



thus account qualitatively correctly for the charging of 
the local level, but the kinks appearing in the I(fJ.) curve 
are due to the incorrectly predicted symmetry breaking 
and are thus again artifacts of the MF solution. 

In the inset of Fig. [5] we have plotted the polarization 
of the current, Pj = (7| — Jj_)/(Jf + JjJ, as a function 
of fi. In th e symmetric case the current is not po larized 




FIG. 5: The current as a function of the bias voltage for 
Ed/U — —0.5 (solid and dashed line) and Ed/U = 0.25 (dot- 
ted line) for T/U = 0.02. Hysteresis can be observed in the 
current in the local moment regime. Inset: the polarization of 
the current for £ d /U = -0.35 and T/U = 0.02. 



due to the electron- hole symmetry, but for small electron- 
hole asymmetries, a hysteresis appears in the polarization 
too. In general, the polarization is finite whenever a local 
moment appears on the impurity and there is no electron- 
hole symmetry. 

To conclude, we have studied the Anderson model out 
of equilibrium in the framework of the scattering for- 
malism combined with a mean-field approximation. This 
method, frequently used in molecular transport calcu- 
lations, incorrectly predicts a magnetic phase transition 
as well as a bias-induced magnetic moment formation, 
accompanied by hysteresis in various physical quantities 
and the coexistence of multiple solutions. The MF ap- 
proach thus fails whenever correlations become impor- 
tant. These artifacts of the mean field approach should 
alert physicists who study transport through strongly 
correlated and magnetic molecules, and urge one to use 
more sophisticated methods that avoid spontaneous sym- 
metry breaking and account for dynamical effects. 
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